clear all
*cap log close
set more off
set seed 603

// macros for figure construction ----------------------
local Q = 2 
local nbin = 15

// which files are we reading in ------------
local file1 = "history_q2"
local file2 = "history_b5"
// -----------------------------------------------------


// ------------------------------------------------------------------------------------------
// first plot: extrapolation estimates -------------

// quadratic tail estimates -------
import delimited using "${est}/`file1'.csv", clear
replace par = trim(upper(par))
keep if par == "Y0"|par=="Y1"
gen z = (par=="Y1")

rename est y 
gen G = g-1
keep G z y ub lb 
gen flag = 1
tempfile est1
save    `est1'

/*
// binscatter tail estimates ---------
import delimited using "${est}/`file2'.csv", clear
replace par = trim(upper(par))
keep if par == "Y0"|par=="Y1"

gen z = .
replace z = 0.025 if par == "Y0"
replace z = 0.975 if par == "Y1"

rename est y 
gen G = g-1
keep G z y ub lb 
tempfile est2 
save    `est2'\
*/

// store quadratic fits ----------------
use "${data}/out/4-main", clear
gen Y =  cite_ny1 
gen D = harsh 
gen jid = officerid 
gen W = covbin1 
gen G = cite_py2

keep  Y D Z W G totfe jid 
order Y D Z W G totfe jid  

// quadratic fit: groups ---------
forval g = 0/1 {
	forval q = 1/`Q' {
		gen g`g'_Z`q' = (G==`g')*(Z^`q')
	}
}
reghdfe Y g*, absorb(fehat=totfe ghat=G) nocons
forval g = 0/1 {
	forval q = 1/`Q' {
		local g`g'_b`q' = _b[g`g'_Z`q']
	}
}

gen sumfe = fehat + ghat 
forval g = 0/1 {
	qui summ sumfe if G == `g'
	local g`g'_b0 = r(mean)
}


// binscatter by groups -------------
binsreg Y Z, absorb(totfe G) by(G) binspos(es) nbins(`nbin') samebinsby savedata(bs) replace

use bs, clear 
gen z = dots_x 
gen y = dots_fit 
qui append using `est1'
*qui append using `est2'
rm bs.dta 

// superimpose quadratic fit --------
gen fit = . 
forval g = 0/1 {
	replace fit = `g`g'_b0' if G == `g'
	forval q = 1/`Q' {
		replace fit = fit + `g`g'_b`q''*(z^`q') if G == `g'
	}
}

// build plot ------------------------
sort G z
#delimit ;
scatter y z if G==0, msymbol(Oh) mcolor(dkgreen%90) msize(medium) yaxis(2) || 
line fit z if G==0, lcolor(dkgreen) lpattern(dash) yaxis(2) ||
rspike ub lb z if G==0 & flag==1, lcolor(dkgreen) yaxis(2) ||
scatter y z if G==0 & flag==1, msymbol(O) mcolor(dkgreen) yaxis(2) ||
scatter y z if G==1, msymbol(Sh) mcolor(purple%90) msize(medium) || 
line fit z if G==1, lcolor(purple) lpattern(dash) ||
rspike ub lb z if G==1 & flag==1, lcolor(purple) ||
scatter y z if G==1 & flag==1, msymbol(S) mcolor(purple) 
graphregion(color(white)) plotregion(lcolor(black) lwidth(medthin)) 
ylab(0.2(0.025)0.3,nogrid axis(2)) ylab(0.4(0.025)0.5,nogrid axis(1)) xlab(,nogrid)
xtitle("Officer Stringency") ytitle("Pr(Reoffend)",axis(2))
legend(ring(0) cols(1) pos(7) order(1 5) 
lab(1 "No Past Offense (Left Axis)") lab(5 "Any Past Offense (Right Axis)")
region(lstyle(border)))  ;
#delimit cr
graph export "${out}/main/history_extrap.pdf", replace
// ------------------------------------------------------------------------------------------




// ------------------------------------------------------------------------------------------
// second plot: levels and gains ---------

// For graphing: groups ------
forval g = 1/2 {
	
	import delimited "${est}/`file1'.csv", clear
	replace par = trim(upper(par))
	keep if g == `g'

	* populate levels and gains ----------------
	do "${code}/utility/create_levels_gains.do"

	drop if mi(y0)
	keep ind y0 y0u y0l y1 y1u y1l delta deltau deltal
	replace ind = 0.2 if _n == 1
	replace ind = 0.5  if _n == 2 
	replace ind = 0.8 if _n == 3

	gen G = `g'-1
	tempfile temp`g'
	save    `temp`g''
	
}



// Compile data -----------
use `temp1', clear
qui append using `temp2'
*replace ind = ind-0.02 if G==3
replace ind = ind-0.02 if G==0


// With axis adjustments ------
#delimit ; 
twoway 
(rbar y1u y1l ind if G==0, color(magenta%20) barwidth(0.02) yaxis(2))
(scatter y0 ind if G==0, msize(medium) msymbol(Dh) mcolor(dkgreen) yaxis(2))
(pcarrow y0 ind y1 ind if G==0, lcolor(magenta) lwidth(medthick) mcolor(magenta) barbsize(medlarge) yaxis(2))
(rbar y1u y1l ind if G==1, color(dkorange%20) barwidth(0.02))
(scatter y0 ind if G==1, msize(medium) msymbol(Sh) mcolor(purple))
(pcarrow y0 ind y1 ind if G==1, lcolor(dkorange) lwidth(medthick) mcolor(dkorange) barbsize(medlarge)),
xscale(r(0 1)) xlabel(0.2 "All Motorists" 0.5 "Treated" 0.8 "Untreated") 
ylab(0.35(0.05)0.5,nogrid axis(1)) ylab(0.15(0.05)0.30,nogrid axis(2)) xlab(,nogrid)
legend(ring(0) pos(7) order(2 5) cols(1) 
label(2 "No Past Offense (Left Axis)") label(5 "Any Past Offense (Right Axis)")
region(lstyle(border))) 
xtitle("") ytitle("",axis(1)) ytitle("Pr(Reoffend)",axis(2)) 
graphregion(color(white)) bgcolor(white) plotregion(lcolor(black) lwidth(medthin)) ;
#delimit cr
graph export "${out}/main/history_est.pdf", replace
// ------------------------------------------------------------------------------------------















